#####################################
# Replication Material for 
# Stefan Müller and Michael Jankowski:
# Do voters really prefer more choice? Determinants of support for personalised electoral systems.
# Journal of Elections, Public Opinion and Parties.
#####################################

## Note: This script creates predicted probability plots based on the Stata output
## from file "02_analyse_hh_hb.do". The Stata file creates csv files with
## predicted probabilities which will be loaded into R and plotted with ggplot2.

# Import function that uses results from Stata do file to create 
# predicted probability plots
source("function_plot.R")

library(tidyverse)
library(srvyr)
library(haven)

# Function für Plots

create_plot <- function(model_name = "pident", 
                        xlab = "", 
                        exclude_do_not_mind = FALSE,
                        trans = FALSE,
                        colors = TRUE){

  dfs <- lapply(1:3, function(x) mimport(paste0(model_name,"_",x))[[1]])
  
  for(i in 1:3) dfs[[i]]$Outcome <- c("Disapproval","Indifference", "Approval")[i]
  
  dfs <- bind_rows(dfs)
  
  dfs$Outcome <- factor(dfs$Outcome, levels = c("Disapproval","Indifference", "Approval"))
  
  if(exclude_do_not_mind == TRUE){
    dfs <- dplyr::filter(dfs, Outcome != "Indifference")
  }
  
  dfs$Predictor <- dfs$Factor1
  
  dfs$Predictor <- sapply(dfs$Predictor, function(x) gsub(" - ","\n-\n",x))
  
  dfs$Predictor[grepl("^Gr.+ne$",dfs$Predictor)] <- "Greens"
  dfs$Predictor[dfs$Predictor == "Waehler\n-\nKein Splitting"] <- "Voter\n-\nno splitting"
  dfs$Predictor[dfs$Predictor == "Waehler\n-\nSplitting"] <- "Voter\n-\nsplitting"
  dfs$Predictor[dfs$Predictor == "Nichtwaehler"] <- "Nonvoter"
  
  if(model_name == "pident" | model_name == "pident_local"){
    dfs$Predictor <- car::recode(dfs$Predictor, "'Linke'='Left'; 'Greens'='Greens';
                                 'Keine'='None';'Andere'='Other'")
    dfs$Predictor <- factor(dfs$Predictor, levels = c("Left", "Greens", "SPD",
                                                      "CDU", "FDP", "None", "Other"))
  }
  
  if(model_name == "age_fine"){
    
    dfs$Predictor <- gsub("age_fine=","", dfs$Predictor)
    dfs$Predictor <- factor(dfs$Predictor, levels = paste(1:11))
    dfs$Predictor <- plyr::revalue(dfs$Predictor, c("1" = "16-17", "2" = "18-20", "3" = "21-24", "4" = "25-29", "5" = "30-34", "6" = "35-39", "7" = "40-44", "8" = "45-49", "9" = "50-59", "10" = "60-69", "11" = "70+"))
  }
  
  if(model_name == "pol_know"){
    dfs$Predictor <- car::recode(dfs$Predictor, "'high level'='High'; 'low level'='Low'")
    dfs$Predictor <- factor(dfs$Predictor, levels = c("Low", "High"))
    
  }
  if(model_name == "pol_int"){
    dfs$Predictor <- car::recode(dfs$Predictor, "'Niedrig'='Low'; 'Mittel'='Medium';'Hoch'='High'")
    
    dfs$Predictor <- factor(dfs$Predictor, levels = c("Low", "Medium", "High"))
  }
  
  if(model_name == "abi"){
    
    dfs$Predictor[dfs$Predictor=="Kein Abitur"] <- "No A levels"
    dfs$Predictor[dfs$Predictor=="Abitur\n-\nkein Studium"] <- "A levels\n-\nno college"
    dfs$Predictor[dfs$Predictor=="Abitur\n-\nStudium"] <- "A levels\n-\ncollege"
    dfs$Predictor <- factor(dfs$Predictor, levels = c("No A levels", "A levels\n-\nno college", "A levels\n-\ncollege"))
  }
  

  if(trans){
    p <- ggplot(dfs, aes(x = Outcome, 
                         y = b, 
                         ymin = llci, 
                         ymax = ulci, 
                         color = Outcome))
  } else{
    p <- ggplot(dfs, aes(x = Predictor, 
                         y = b, 
                         ymin = llci, 
                         ymax = ulci, 
                         color = Predictor))
  }
  
  p <- p + geom_pointrange(size = .75) +
    scale_y_continuous(limits = c(0, max(dfs$ulci)), breaks = c(seq(0, max(dfs$ulci), 0.1))) +
    xlab(xlab) +
    facet_wrap(~ Outcome, ncol = 3) +
    ylab("Predicted Probability") +
    guides(color = FALSE) +
    theme_bw(base_size = 16) +
    theme(
      panel.background = ggplot2::element_blank(),
      panel.grid.major.x = element_blank(),
      panel.grid.minor.x = element_blank(),
      panel.grid.minor.y = element_blank(),
      strip.background = element_rect(colour = "white", fill = "white"),
      axis.text = element_text(colour="black"),
      axis.title = element_text(colour="black"),
      strip.text = element_text(colour="black", size = 16),
      plot.background = element_blank(),
      panel.grid.major.y = element_blank()
    )
  if(colors) p <- p + scale_color_manual(values = colorRampPalette(c('#132145', '#89c7ff'))(length(unique(dfs$Predictor))))
  if(!colors) p <- p + scale_color_manual(values = rep("black", length(unique(dfs$Predictor))))
  p
}



############################
# Summary statistics (Figure 1) ----
############################

df <- haven::read_dta("hh_hb_11_15_prepared.dta")


# Summarise data using srvyr package
df_summary <- df %>% 
    filter(!is.na(pref_law)) %>% 
    srvyr::as_survey_design() %>% 
    group_by(state, year, election_factor, pref_law_factor) %>% 
    summarize(prop = srvyr::survey_mean(na.rm = TRUE)) 

df_summary$pref_law_factor <- factor(df_summary$pref_law_factor,
                                     levels = c("Disapproval", "Indifference", "Approval"))

dodge <- position_dodge(width = 0.7)

ggplot(data = df_summary, 
       aes(x = pref_law_factor, 
           y = prop, 
           fill = year)) +
    geom_bar(stat = "identity",  width = 0.6, position = dodge) +
    labs(x = "Attitude towards the electoral system",
         y = "Percentage of respondents") +
    geom_errorbar(aes(ymin = prop - 2 * prop_se,
                      ymax = prop + 2 * prop_se),
                  position = dodge, width = 0.3) +
    scale_fill_grey(start = 0.7, end = 0.3, name = "Year") + 
    facet_wrap(~ state, scales = "free_x") +
    scale_y_continuous(labels = scales::percent) +
    theme_bw(base_size = 16) +
    theme(
        legend.position = "bottom",
        panel.background = ggplot2::element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        strip.background = element_rect(colour = "white", fill = "white"),
        axis.text = element_text(colour="black"),
        axis.title = element_text(colour="black"),
        strip.text = element_text(colour="black", size = 16),
        plot.background = element_blank(),
        panel.grid.major.y = element_blank()
    )


############################
# Party identification (Figure 2) ----
############################


colors_party <- c("#96276E", "#4BA345","#DB4240", "#373737",  "#F6BB00",  "grey70",   "lightblue")

pident <- create_plot("pident", "Party identification") +
  scale_color_manual(values = colors_party) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) 

pident

############################
# Political interest (Figure 3) ----
############################

pol_int <- create_plot("pol_int", "Political interest")

pol_int


############################
# Age (Figure 4) ----
############################

age <- create_plot("age", "Age (grouped)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) 

age

########################################################
########################################################
### Figures for Supplementary Material
########################################################
########################################################


############################
# Political knowledge (Figure A1) ----
############################

pol_know <- create_plot("pol_know", "Political knowledge")

pol_know




############################
# Satisfaction (Figure A2) ----
############################

dfs <- lapply(1:3, function(x) mimport(paste0("satisfaction_", x))[[1]])
for(i in 1:3) dfs[[i]]$Outcome <- c("Disapproval", "Indifference", "Approval")[i]
dfs <- bind_rows(dfs)

dfs$Outcome <- factor(dfs$Outcome, levels = c("Disapproval","Indifference", "Approval"))

satis <- ggplot(dfs, aes(x = satisfaction, y = b, ymin = llci, ymax = ulci, fill = Outcome)) +
  geom_ribbon(alpha = .5, color = 0) +
  geom_line(color = "black") +
  facet_wrap(~ Outcome, ncol = 3) +
  xlab("Satisfaction with party") +
  ylab("Predicted Probability") +
  scale_x_continuous(breaks = 0:10) +
  scale_color_manual(values = c("#991337","#1d7eed","#5da352")) + 
  scale_fill_manual(values = c("#991337","#1d7eed","#5da352")) +
  theme_bw(base_size = 16) +
  theme(
    legend.position = "none",
    panel.background = ggplot2::element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    strip.background = element_rect(colour = "white", fill = "white"),
    axis.text = element_text(colour="black"),
    axis.title = element_text(colour="black"),
    strip.text = element_text(colour="black"),
    plot.background = element_blank(),
    panel.grid.major.y = element_line(linetype = "dotted"))

satis



############################
# Age Fine Grained (Figure A3) ----
############################


age_fine <- create_plot("age_fine", "Age (grouped)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) 

age_fine 


############################
# Education (Figure A4) ----
############################

education <- create_plot("abi", "Formal education")

education


############################
# Voter type (Figure A5) ----
############################

splitting <- create_plot("splitting", "Voter type")

splitting


############################
# Party identification local (Figure A6) ----
############################

pident_local <- create_plot("pident_local", "Party identification (local level)") +
  scale_color_manual(values = colors_party_modified) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) 

pident_local


######################
# Party ratings (Figure A7) ----
######################

read_mod <- function(model_name, party){
dfs <- lapply(1:3, function(x) mimport(paste0(model_name,"_",x))[[1]])
for(i in 1:3) dfs[[i]]$Outcome <- c("Ablehnung","Egal","Zustimmung")[i]
dfs <- bind_rows(dfs)
dfs$Partei <- party
dfs
}

parties <- c("CDU", "SPD", "left", "green", "FDP")

mods <- plyr::ldply(parties, function(x) read_mod(paste0("rating_",tolower(x)), x))

mods$Partei[mods$Partei == "green"] <- "Greens"
mods$Partei[mods$Partei == "left"] <- "Left"

mods$rate <- as.numeric(apply(mods[,grepl("rate_", colnames(mods))], 1, function(x) x[!is.na(x)]))

mods$Outcome <- car::recode(mods$Outcome, "'Ablehnung'='Disapproval'; 'Zustimmung'='Approval';
                            'Egal'='Indifference'")

mods$Outcome <- factor(mods$Outcome, levels = c("Disapproval", "Indifference", "Approval"))

rating_plot <- ggplot(mods, aes(x = rate, y = b, ymin = llci, ymax = ulci)) +
  geom_ribbon(aes(fill = Partei), color = 0, alpha = 0.5) + 
  geom_line(color = "black") +
  facet_grid(Partei ~ Outcome) +
  ylab("Predicted probability") + 
  xlab("Approval ranking of party") +
  theme_bw(base_size = 16) +
  theme(
    panel.background = ggplot2::element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    strip.background = element_rect(colour = "white", fill = "white"),
    axis.text = element_text(colour="black"),
    axis.title = element_text(colour="black"),
    strip.text = element_text(colour="black"),
    plot.background = element_blank(),
    panel.grid.major.y = element_line(linetype = "dotted")
  ) +
  scale_color_manual(values = colors_party_modified[c(4,5,2,1,3)], guide = FALSE) +
  scale_fill_manual(values = colors_party_modified[c(4,5,2,1,3)], guide = FALSE) +
  scale_x_continuous(breaks = 0:10)

rating_plot



############################
# Satisfaction interaction (Figure A8) ----
############################

dfs <- lapply(1:3, function(x) mimport(paste0("supporters_int_", x))[[1]])
for(i in 1:3) dfs[[i]]$Outcome <- c("Disapproval", "Indifference", "Approval")[i]
dfs <- bind_rows(dfs)

dfs$Outcome <- factor(dfs$Outcome, levels = c("Disapproval","Indifference", "Approval"))

dfs$supporter <- ifelse(dfs$Factor1 == "supporters=1", "Left Party/Greens", "Other Parties")

satis_int <- ggplot(dfs, aes(x = satisfaction, y = b, ymin = llci, ymax = ulci, fill = supporter)) +
  geom_ribbon(alpha = .5, color = 0) +
  geom_line(color = "black") +
  facet_wrap(~ Outcome, ncol = 3) +
  xlab("Satisfaction with party") +
  ylab("Predicted Probability") +
  scale_x_continuous(breaks = 3:10) +
  theme_bw(base_size = 16) +
  theme(
    legend.position = "bottom",
    panel.background = ggplot2::element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    strip.background = element_rect(colour = "white", fill = "white"),
    axis.text = element_text(colour="black"),
    axis.title = element_text(colour="black"),
    strip.text = element_text(colour="black"),
    plot.background = element_blank(),
    panel.grid.major.y = element_line(linetype = "dotted")) +
  guides(fill = guide_legend(title = "Party"))

satis_int



create_plot_sep <- function(model_name = "pident", 
                        xlab = "", 
                        exclude_do_not_mind = FALSE,
                        trans = FALSE,
                        colors = TRUE){
  
  dfs_1 <- lapply(1:3, function(x) mimport(paste0(model_name,"_hb_2011_",x))[[1]])
  dfs_2 <- lapply(1:3, function(x) mimport(paste0(model_name,"_hb_2015_",x))[[1]])
  dfs_3 <- lapply(1:3, function(x) mimport(paste0(model_name,"_hh_2011_",x))[[1]])
  dfs_4 <- lapply(1:3, function(x) mimport(paste0(model_name,"_hh_2015_",x))[[1]])
  
  for(i in 1:3) dfs_1[[i]]$Outcome <- c("Disapproval","Indifference", "Approval")[i]
  for(i in 1:3) dfs_2[[i]]$Outcome <- c("Disapproval","Indifference", "Approval")[i]
  for(i in 1:3) dfs_3[[i]]$Outcome <- c("Disapproval","Indifference", "Approval")[i]
  for(i in 1:3) dfs_4[[i]]$Outcome <- c("Disapproval","Indifference", "Approval")[i]
  
  dfs_1 <- bind_rows(dfs_1)
  dfs_1$election <- "Bremen 2011"
  dfs_2 <- bind_rows(dfs_2)
  dfs_2$election <- "Bremen 2015"
  dfs_3 <- bind_rows(dfs_3)
  dfs_3$election <- "Hamburg 2011"
  dfs_4 <- bind_rows(dfs_4)
  dfs_4$election <- "Hamburg 2015"
  
  dfs <- bind_rows(list(dfs_1, dfs_2, dfs_3, dfs_4))
  
  dfs$Outcome <- factor(dfs$Outcome, levels = c("Disapproval","Indifference", "Approval"))
  
  if(exclude_do_not_mind == TRUE){
    dfs <- dplyr::filter(dfs, Outcome != "Indifference")
  }
  
  dfs$Predictor <- dfs$Factor1
  
  dfs$Predictor <- sapply(dfs$Predictor, function(x) gsub(" - ","\n-\n",x))
  
  #dfs <- subset(dfs, Predictor != "Andere")
  
  dfs$Predictor[grepl("^Gr.+ne$",dfs$Predictor)] <- "Greens"
  dfs$Predictor[dfs$Predictor == "Waehler\n-\nKein Splitting"] <- "Voter\n-\nno splitting"
  dfs$Predictor[dfs$Predictor == "Waehler\n-\nSplitting"] <- "Voter\n-\nsplitting"
  dfs$Predictor[dfs$Predictor == "Nichtwaehler"] <- "Nonvoter"
  
  if(model_name == "pident" | model_name == "pident_local"){
    dfs$Predictor <- car::recode(dfs$Predictor, "'Linke'='Left'; 'Greens'='Greens';
                                 'Keine'='None';'Andere'='Other'")
    dfs$Predictor <- factor(dfs$Predictor, levels = c("Left", "Greens", "SPD",
                                                      "CDU", "FDP", "None", "Other"))
    dfs <- subset(dfs, Predictor != "Other")
  }
  
  if(model_name == "age_fine"){
    
    dfs$Predictor <- gsub("age_fine=","", dfs$Predictor)
    dfs$Predictor <- factor(dfs$Predictor, levels = paste(1:11))
    dfs$Predictor <- plyr::revalue(dfs$Predictor, c("1" = "16-17", "2" = "18-20", "3" = "21-24", "4" = "25-29", "5" = "30-34", "6" = "35-39", "7" = "40-44", "8" = "45-49", "9" = "50-59", "10" = "60-69", "11" = "70+"))
  }
  
  if(model_name == "pol_know"){
    dfs$Predictor <- car::recode(dfs$Predictor, "'high level'='High'; 'low level'='Low'")
    dfs$Predictor <- factor(dfs$Predictor, levels = c("Low", "High"))
    
  }
  if(model_name == "pol_int"){
    dfs$Predictor <- car::recode(dfs$Predictor, "'Niedrig'='Low'; 'Mittel'='Medium';'Hoch'='High'")
    
    dfs$Predictor <- factor(dfs$Predictor, levels = c("Low", "Medium", "High"))
  }
  
  if(model_name == "abi"){
    
    dfs$Predictor[dfs$Predictor=="Kein Abitur"] <- "No A levels"
    dfs$Predictor[dfs$Predictor=="Abitur\n-\nkein Studium"] <- "A levels\n-\nno college"
    dfs$Predictor[dfs$Predictor=="Abitur\n-\nStudium"] <- "A levels\n-\ncollege"
    dfs$Predictor <- factor(dfs$Predictor, levels = c("No A levels", "A levels\n-\nno college", "A levels\n-\ncollege"))
  }
  
  
  if(trans){
    p <- ggplot(dfs, aes(x = Outcome, 
                         y = b, 
                         ymin = llci, 
                         ymax = ulci, 
                         color = Outcome))
  } else{
    p <- ggplot(dfs, aes(x = Predictor, 
                         y = b, 
                         ymin = llci, 
                         ymax = ulci, 
                         color = Predictor,
                         group = election,
                         shape = election))
  }
  
  p <- p + geom_pointrange(size = .75, position = position_dodge(width = 0.7)) +
    scale_y_continuous(limits = c(0, max(dfs$ulci)), breaks = c(seq(0, max(dfs$ulci), 0.1))) +
    xlab(xlab) +
    facet_wrap(~ Outcome, ncol = 3) +
    ylab("Predicted Probability") +
    guides(color = FALSE, shape = guide_legend(title = "Election")) +
    theme_bw(base_size = 24) +
    theme(
      panel.background = ggplot2::element_blank(),
      panel.grid.major.x = element_blank(),
      panel.grid.minor.x = element_blank(),
      panel.grid.minor.y = element_blank(),
      strip.background = element_rect(colour = "white", fill = "white"),
      axis.text = element_text(colour="black"),
      axis.title = element_text(colour="black"),
      strip.text = element_text(colour="black", size = 28),
      plot.background = element_blank(),
      panel.grid.major.y = element_blank()
    ) +
    scale_shape_manual(values = c(0,1,15,16)) +
    theme(legend.position = "bottom")
  if(colors) p <- p + scale_color_manual(values = colorRampPalette(c('#132145', '#89c7ff'))(length(unique(dfs$Predictor))))
  if(!colors) p <- p + scale_color_manual(values = rep("black", length(unique(dfs$Predictor))))
  p
}



############################
# Party identificaiton (separate) (Figure A9) ----
############################

pident_sep <- create_plot_sep("pident", "Party identification") +
  scale_color_manual(values = colors_party) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) 

pident_sep 

###############################
# Political interest (separate) (Figure A10) ----
###############################

pol_int_sep <- create_plot_sep("pol_int", "Political interest")

pol_int_sep


############################
# Age (separate) (Figure A11) ----
############################

age_sep <- create_plot_sep("age", "Age (grouped)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) 

age_sep


############################
# Education (separate) (Figure A12) --
############################

education_sep <- create_plot_sep("abi", "Formal education")

education_sep

##################################
# Voter type (separate) (Figure A13) ----
##################################

splitting_sep <- create_plot_sep("splitting", "Voter type")

splitting_sep

